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ABSTRACT 

We present spherical, non-rotating, isotropic models of early-type galaxies with stellar and dark- 
matter components both described by deprojected Sersic density profiles, and prove that they represent 
physically admissible stable systems. Using empirical correlations and recent results of N-body sim- 
ulations, all the free parameters of the models are expressed as functions of one single quantity: the 
total (B-band) luminosity of the stellar component. 

We analyze how to perform discrete N-body realizations of Sersic models. To this end, an optimal 
smoothing length is derived, defined as the softening parameter minimizing the error on the gravita- 
tional potential for the deprojected Sersic model. It is shown to depend on the Sersic index n and on 
the number of particles of the N-body realization. 

A software code allowing the computations of the relevant quantities of one- and two-component 
Sersic models is provided. Both the code and the results of the present work are primarily intended 
as tools to perform N-body simulations of early-type galaxies, where the structural non-homology of 
these systems (i.e. the variation of the shape parameter along the galaxy sequence) might be taken 
into account. 

Subject headings: Galaxies - Astronomical Techniques 



1. INTRODUCTION 



Merging of red-sequence 
important channel for the 
early- type galaxies (ETGs). 
ers have been observed to 
the 



galaxies might be an 
formation of massive 
Such dry merg- 
take place and have 
an impact on the popula tion of ETGs at both 
low (u p to z ^ 0.3; IWhitaker fc van DokkumI 
120081 iMasjedi. Hogg fc Blanto^ I2008D . and in- 
termedi ate redshift, in cluster and field environ- 
ments (Ivan Do kkum et al.' '19991 Ivan DokkumI I2OO5I : 
iTran et all 1200 5: BcU ct al. 200^. A further evidence 
comes from the fact that the stellar mass on the red 
sequenc e has been found to be nearly d oubled from 
z ~ 1 (jZucca et al.l 120061 : iBell et al.l l2004f ) on, imply- 
ing that at least some red galaxies must be formed 
from mer ging systerns tha t are either very dusty or 
gas-poor (|Faber et al.ll2005l ). K-band selected samples 
also revealed a substantial population of old, passively 
evolving, massive ETGs already in place at 1 < z < 2, 
with luminosity and stellar mass function s evolving 
only weakly up to z ^ 0.8 — 1 jC imatti et al.l f200l : 
iBundv et al]l2006l : ICimatti. Daddi fc Renzinill2006f ). 
From the theoretical viewpoint, dry mergers are also 

expecte d to play a major ro le. Using semi-analytical 

models, iKhochfar fc Burker^ (|2003l ) found that a large 
fraction of present-day ETGs are indeed formed by 
merging bulge-dominated systems and that the fraction 
of spheroidal mergers increases with luminosity, with 
massive ETGs being formed by nearly d issipationless 
events. As shown bv iDe Lucia et al] (|2006D . more mas- 
sive ETGs are expected to be built up of several stellar 
pieces, with the number of effective stellar progenitors 
increasing up to five for the most massive galaxies. 
On the other hand, hydro-dynamical simulations have 



also shown that accretion of smaller disk-dominated 
galaxies (in the mass ratio of 1:10) could also have 
an important role in the evolution of massive ETGs, 
explaining the pres ence of t he tidal debris observed at 
z - (Feldmann. Maver fc CarollQ .2008). 

To constrain the role of dry mergers in galaxy 
formation, it is of importance to perform merging 
simulations of spheroidal systems, comparing the 
properties of merger remnants to observations. So 
far, merging simulations of ETGs have been mostly 
used to constrain the origin of the empirical corre- 
lations among gala xy's obse rved quantities, such as 
the Faber-Jackson (iFaber fc J ackson 1976L hereafter 
FJ) , the Kormendy (jKormendvl Il977l hereafter KR), 
and the Fundamental Plane (|Diorgovski fc DavisI 
Il987l hereafter FP) relations. The impact of 
dry merging has been investigated in several 
works ( e.g. ICapelatq. de Carvalho fc Carlberd 

"2004; 



'1995'; 'Dantas et al. '2003: 'E vstigneeva et a l' ' 



Nipoti, Londrillo fc Ciotti 2003). They have all agreed 
that dissipationless merging is able to move galax- 
ies along the FP. But it is not clear if dry mergers 
are al so able to preserve other observed correla- 
tions (Bovlan-Kolchin, Chung-Pci fc Eliot' '2006). For 
instance, Nipoti, Londrillo fc Ciotti (2003) found that 
the products of repeated merging of gas-free galaxies 
are characterized by an unrealistically large effective 
radius and a mass-ind e pende nt velocity dispersion, 
while lEvstigneeva et ahl (|2004f ) found that only the 
merging of massive galaxies, that lie on the KR, leads 
to end-products that still follow that relation. 
In previous works, merging simulations have been 
performed by means of ETG's models where the stellar 
component is described by simple analytic density 
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laws, such as the King or the Hernquist profiles. This 
approach implicitly neglects one key observational 
feature: the structural non-homo logy of the ETG pop- 
ulation (jGraham fc Collesill997D . It is well established 
that the observed light profiles of ETGs deviate from 
a pur e r^/** law, being better described by the Sersic 
([1961) model (ICaon. CapaccioH fc D'Onofriol [1991 



D'Onofrio . CapaccioH fc Cam] [1991 \g raham et alj 
1996). The Sersic index (shape parameter), n, measuring 
the steepness of the light profile, changes systematically 
along the galaxy sequence, the more luminous galaxies 
having higher n. Moreover, the shape parameter also 
correlates with other observed properties of ETGs, such 
as the effe ctive paramet ers and the central velocity 
dispersion (jGraham! 1200 ^. as expected in view of the 
correlation of n with the luminosity. Different values of 
n correspond to physical systems that differ significantly 
in their phase-space density structure, with higher 
Sersic indices describing galaxies whose light profile is 
significantly more concentrated toward the center, with 
an extended low surface brightness halo. Thus, merging 
systems with different n's might lead to a different evo- 
lution of the phase-space density of merging remnants 
with respect to that of "homologous" King/Hernquist 
models. For what concerns dark matter haloes, previous 
simulations have usually ado pted either the Navarro- 
Frenk -White (NEW) pro file (iNavarro. Erenk fc White! 
!l995l ) or th e Hernquist (1990!) p rofile. However, as 
shown by !Merritt et al.! (j2005! ) and !Merritt et all 
()2006l ) (hereafter MGM06), galaxy- and cluster-sized 
halos are actuall y better desc ribed by using either the 
Einasto's model (lEinastol!l968D or the Prugniel fc Simien 
model (!Prugniel fc Simien 19971) rather than a NFW-like 
profile (jNavarro. Erenk fc White! [19951 ). The Einasto's 
model is identical in functional form to the Sersic model, 
but is used to describe the deprojected (rather than the 
projected) density profile, while the Prugniel fc Simien 
model is an an alytic approximation to the deprojected 
Sersic profile. !Merritt et all (|2005D (hereafter MNL05) 
and MGM06 found that the deprojected Sersic model 
(i.e. the Prugniel & Simien model) provides a better 
fit to the projected mass density profile of simulated 
dark-matter halos, with a Sersic index value of n ~ 3 for 
galaxy-sized dark-matter halos. 

Hence, the deprojected Sersic model seems able to de- 
scribe both the stellar and dark matter components of 
ETGs. Driven by that, we present here new simple 
models of ETGs, where both components follow the de- 
projected Sersic law. Hereafter, we refer to these mod- 
els as double Sersic (5'^) models. The models describe 
spherical, non-rotating, isotropic systems, and are in- 
tended as a tool to perform N-body simulations of ETGs. 
In a companion contribution (Coppola et al. 2009b, in 
preparation), we use the models to investigate how 
dissipation-less (major and minor) mergers affect the 
structural properties of ETGs, such as the shape of their 
light profile and their stellar population gradients. The 
present paper aims at: (i) describing the main charac- 
teristics of the S"^ models, by deriving the corresponding 
potential-density pair and distribution function (Sec. 2), 
and discussing their physical consistency and stability 
(Sec. 3); (ii) describing how to perform discrete N-body 
realization of the models, by adopting an optimal gravi- 



tational smoothing length for simulation codes (Sec. 4); 
(iii) giving a set of recipes to fix all the free model param- 
eters (Sec. 5); (iv) providing the software code to com- 
pute dynamical/structural properties of both the one- 
and two-component Sersic models. Summary and dis- 
cussion are drawn in Sec. 6. 

2. THE DOUBLE SERSIC (S^) MODEL 

2.1. The deprojected Sersic model 

The surface brightness profile of ETGs, 
I{R), is accurately described by the Ser- 
sic law (!Capaccioli. Caon fc D'Onofrio! 

1991 !Caon. Capaccio H fc D'Onofriol 



1991 



D'Onofrio. Capaccioh fc Caon 1994): 



I{R;n) = /oexp -b(R/R. 



^l/n 



(1) 



where Jo is the central surface brightness, R is the (equiv- 
alent) projected distance to the galaxy center, n is the 
Sersic index (shape parameter), and & is a function of n, 
defined in such a way that Rp , is the eff ective (half-light) 
radius of the galaxy (|Ciottil!l99ll !l999! ). The quantity b 

is approximated at better than 1% by the relation b ^ 

exp [0 .6950 + ln(n) - 0.1789/n] ([Lima Neto. Gerbal fc Marqued 
[l999l) . 

For a spherical system, under the assumption that the 
stellar mass-to-light ratio, M^^/L, does not change with 
radius, the spatial mass density profile of the stellar com- 
pone nt, p, , is obtained by solvin g the Abel integral equa- 
tion (jBinnev fc Tremaine! [19881 ) . 



Pl (r) = - 



TT L 



dl 



dR ViF 



dR 



(2) 



where r is the distance to the galaxy center. Setting 
u — r/R and inserting Eq. [l]into the Abel equation, one 
obtains the following expression: 



exp[ 



-6a;i/"w-i/"l du 



where x = r/Re^ is the distance to the galaxy center 
in units of i?ei , p(a;; n) is the dimensionless deprojected 
density profile, and po^ = M^/Rl^ ■ 62"/(27m r(2n)) is 
the scaling factor of the stellar density profile. Here, F 
denotes the complete gamma function, a nd the expres- 
sion of poi, is obtained by using eq. 4 of [Ciottil (|1999D . 
which gives the total luminosity of the Sersic model as a 
function of /q, Re^i ^-i^d From Eq. [3l one obtains the 
mass profile: 



M^{r-n)^MQ^M{x-n) = 



Oil 



62" 7o (l_y2)l/2 
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2n-f l,5f- 



l/n 



du, (4) 



where M is the dimensionless mass profile, and Afoi, = 
A462"/(27rn F(2n)) is the scaling factor of (r). From 
the Laplace equation, one finds the following expression 
for the gravitational potential: 



M (x; n) 



'^^-^ «(i-?/^)-*7^« + i,&(9 j*) 
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where ip{x; n) is the dimensionless gravitational poten- 
tial, and ipoL = GM^/Re^ is the corresponding scaling 
factor. As shown in Sees. 2.2 and 2.3, the above equa- 
tions provide the essential ingredients to construct the 
S*^ models. 

We notice that, due to the existence of radial gradients 
in stellar populati on properties (such as age and metallic- 
ity) of ETGs (e.g. lPeletier et al.lll990D . the assumption of 
a constant mass-to-light ratio {M^/L{r) = const.) might 
not actually reflect the physical properties of early-type 
systems. As discussed in Sec. [H considering the observa- 
tional results on age and metallicity gradients in ETGs, 
the M^/L is expected to vary significantly with galaxy 
radius (up to ~ 50%) at optical wavebands (B-band). 
However, the variation is significantly reduced, becoming 
consistent with zero within observational uncertainties, 
at Near-Infrared (NIR) wavebands. According to that, 
we implicitly assume here that the parameters Re^ and 
n, entering the normalization factors of the potential- 
density pair of the S*^ models, are those describing the 
NIR profile of ETGs. In Sec.[31 we describe how to derive 
the free parameters of the 5^ models according to this 
assumption. 

The deprojection of the Sersic law has b een al- 

(ICiottil 

d 



and mass profile are then obtained from the equations: 



ready presented 



works 



in several 

Prugniel fc SimienI |1997| : iMazure fc Capelat^ 



1991 
200a 



Terzic & Graham 2005). Following iMellier fc Mathied 
(I1987D . IP lugniel fc Simieni ()1997l ) provided an analyt- 



ical approximation to the spatial density profile of the 
i?^/" model (Eq. [S]). iLima Neto. Gerbal fc Marguezl 
()1999l ) showed that the Prugniel fc Simien approxima- 
tion reproduces the deprojected Sersic profile with an 
accuracy better than 5%, in the radial range of 10~^ to 
10'^i?ei,, for Sersic indices between n ^ 0.5 and n ~ 10. 
Th e Prugniel fc Sim i en mo del has been also adopted 
by iTerzic fc Grahaiiil (|2005[ ) to present one-component 
Sersic models of ETGs with power-law cores. Exact 
solutions to the deprojection of the R}^"" m odel have 
been provided by IMazure fc Capelatol ()2002l), in terms 
of the so-called Meijer G functions, while iCiottil ()1991f ) 
presented exact numerical expressions for the mass, 
gravitational potential, and central velocity dispersion 
of the one-component Sersic model. In the present work, 
we report a concise reference to the integral equations 
that define the density-potential pair, the mass profile 
and the distribution function of the deprojected Sersic 
law. All the quantities characterizing the Sersic model 
can be numerically computed by using a set of publicly 
available Fortran programs (see App. A). 



2.2. The dark matter Sersic model 

MNL05 and MGM06 found that the deprojected Ser- 
sic law provides a better fit to the density profile of dark 
matter halos than the NEW law. MNL05 found that a 
Sersic index value of n = 3.00 ±0.17 is required to fit the 
profile of galaxy-sized halos. On the other hand, MGM06 
fitted the Prugniel fc Simien model to the density pro- 
files of galaxy-sized halos, finding a best-fitting value of 
n ^ 3.59 ± 0.65. Considering the lower uncertainty of 
the MNL05 estimate, we describe the dark matter com- 
ponent of the models with a deprojected Sersic model 
having n — i. The corresponding density-potential pair 



M„{r)=tiMo^ Af — ;n = 3 

, ^ fJ. ^ ( X 

(Po{r) = — (POl ^[ — ;n = 3 



(6) 
(7) 
(8) 



where the dimensionless density-potential pair (p, (p) and 

the dimensionless mass profile M are obtained by setting 
n = 3 in Eqs. [3l |4] and [5l respectively. Here, we have 
denoted as = M^/M^ the ratio of the total halo mass, 
M^, to the total stellar mass M^, and a;^ = Re^/Rei, 
the ratio of the (projected) effective radii of the dark 
matter and stellar components. 

We notice that although we fix here the shape param- 
eter value of the dark-matter component, the models 
could be directly generalized to the case where the Ser- 
sic index of the halo component changes with its mass ^ . 
Such a dependece is somewhat suggested by the results 
of MNL05 and MGM06, who found that cluster-sized 
halos [Mjj ^ IQ^^Mq) are better described with Ser- 
sic index values of 2.38 ± 0.25 and - 2.89 ± 0.49, re- 
spectively, these values being systematically smaller than 
those obtained for galaxy-sized halos. However, one 
should notice that, when fitting dwarf-sized dark mat- 
ter halos (M„ - 101°Mq), MNL05 found a best-fitting 
Sersic index value of 3.11 ± 0.05, which is fully consis- 
tent with that of 3.00 ± 0.17 found for galaxy-sized halos 
[M^ ~ IQ^'^Mq). Hence, current results seem to suggest 
a very similar Sersic index value of ~ 3 for galaxy-sized 
halos of different masses, supporting our assumption of 
a fixed n value. 

2.3. Density-potential pair and Distribution Function 

The total mass density profile is obtained by adding up 
the profiles of the stellar and dark matter components: 



pW = Pl+Pd ^ PQl 



p{x-n) + A- p{ — ;3 

xt \ x„ 



(9) 



From the linearity of the Laplace equation, the total 
gravitational potential is equal to ip{r) = ipj^+(p^ , where 
(p^ and (p^ are obtained from Eqs. [5] and [H A simi- 
lar expression can also be obtained for the mass profile, 
combining Eqs. [4] and [Tj We note that the global density- 
potential pair and the mass profile are completely defined 
from five parameters, which are the dimensional quanti- 
ties and Re^ , and the dimension-less parameters x^ , 
/i, and n. 

The distribution function of a stationary, spherical, 
isotropic system depends only on the binding energy 
E and is uniquely defined by the d ensity-potential pair 
through the Eddington formula ^ (jBinnev fc Tremaind 

^ To this aim, one should change Eqs.[6l[7l and[8l by replacing 
the value of n = 3 with a different Sersic index of the dark matter 
halo , and derive the distribution function of the model accordingly 
Sec. ED 

^ As usually done, we write the Eddington formula by adopting 



natural units, where Mj^ 
the gravitational constant. 



1, Re 



1, and G = 1, with G being 
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[Ml): 



71^ [Jo d^T^^ ^ 7S [d^/ *=oJ 

(10) 

where ^(r) = —</'('') + ^^'^ ^ = + fo is the rela- 
tive binding energy, with tpn being a suitably defined con- 
stant (see Binney fc Tremainc 1988). For the S'^ models, 
the global potential and density profiles are proportional 
to the dimensional factors (/?0i, (see Eqs. 5 and 8) and po^ 
(see Eqs. 3 and 6). Hence, using Eq. Bl in App. [B] one 
finds that, unless of a scaling factor depending on 



necessary condition for ^ > is given by: 
g{r]n,p,x„) = 



■d'p 


fd^\ 


dp 






dr dr"^ 



>0. (11) 



For the two-component Sersic models, g{r) is derived nu- 
merically as described in App. [Bj Fig. [5] plots g{r) as a 
function of r for the same sets of n, p, and values 
as in Fig. [1] The condition g{r) > is always verified, 
proving the stability of 5*^ models. 

4. PHYSICAL SCALES 

There are five free parameters that completely charac- 
terize the S'^ model, i.e. the mass of the stellar compo- 
nent, Mj^ , its effective radius, R^j^ , the Sersic index of the 



and Re^ , the J{£) is determined by the three dimension- ^^^^la^ component, n, the mass of the dark matter halo, 
less parameter s x,^, / ' ' - r^__ i_ 
Sersic models ([Ciotti 



p, and 
j fl99l 



n. As for the case of single 
one can show that the sec- 
ond term on the right side of Eq. [10] is always equal to 
zero for all possible values of x^-,, p, and n. In fact, one 
can write {dp/ d'ii)^,^^ = \\mr^oo{dp/dr){dr/d'^). For 
r — > oo, the first derivative of the gravitational potential 
decreases as r~^, while the first deriva tive of the d ensity 
decreases exponentially (see eq. 8 of ICiottil |1991D . im- 
plying that \m\r^oo{dp/dr){dr/d^) — 0. In App. [B] we 
report in detail how to calculate the distribution func- 

,2 

tion by expressing the function in terms of the first 
and second derivatives of p, the gravitational potential 
1^, and the mass profile of the dark matter and stellar 
components. 



3. PHYSICAL CONSISTENCY AND STABILITY 

The Eddington inversion does not guarantee that the 
distribution function is a physically admissible stationary 
solution of the Boltzmann equation. To this effect, for a 
given density-potential pair, one has to show that /(£) is 
non-negative for all positi ve values of the relative bind- 
ing energy. As shown bv ICiottil ([1991i) . one-component 
spherical, non-rotating, isotropic Sersic models are al- 
ways physically admissible, while in the anisotropic case, 
a minimum anisotropy radius exists for the model to be 
admiss ible, with this radius d epending on the Sersic in- 
dex n (jCiotti fc Lanzonil[T997[ ). 

The distribution function of the S"^ models is computed 
by numerical integration of the Eddington formula, as 
described in App. [BJ Fig. [1] plots the f{£) for different 
values of the free parameters n, p and . The value of 
p is varied in the range of zero - no dark matter halo - 
to a value of 10^, where the stellar component is negligi- 
ble and the system is completely dark matter dominated. 
We consider values of x^ from 0.1 to 10^, corresponding 
to the two extreme cases where the dark matter com- 
ponent is either more concentrated or significantly more 
extended than the luminous one. For all combinations 
of and p, different values of n are plotted. We find 
that for positive values of the relative binding energy the 
condition f{£) > is always fulfilled, implying that the 
models are physically admissible. 

To analyze the stabili t y of t he two-component Sersic 
models, following ICiottil ([19911 ). we study the sign of the 
first derivative of the dis tribution function. Acco rding 
to Antonovs theorem fsee lBinnev fc Tremainel[i98^ pag. 
306), if ^ > 0, the system is stable against both radial 
and non-radial perturbations. As shown in App. [B] a 



M^, and the corresponding effective radius, Ren- Al- 
ternatively, one can use the dimensional quantities, 
and Rf,^ , and the dimension-less parameters x^ , p, and n 



defined in Sec. 12.21 Here, we describe some recipes to ex- 
press all the free parameters as a function of one single 
quantity, the absolute luminosity of the stellar compo- 
nent. This procedure is intended as an handy tool to use 
the iS*^ models in merging simulations of ETGs. We refer 
to absolute magnitudes in the B band, Mb, since most 
of the relations we use in the following are expressed in 
that band. In the following, magnitudes are expressed 
with respect to the Vega system. 

The quantity R,,^ is related to the total lumi- 
nosity by the Kormendy relation ([Kormendvi 119771 : 
ICapaccioli. Caon fc D'0^^ofriol[T99l : 



log Re 



(12) 



where Re^ B is the galaxy effective radius in the B-band 
and is the mean effective surface brigthness inside 
Re^ g- Expressing Re^ g units of kpc, one has 

^ -5 \og{Re ^ ) - Mb + 25 + 2.5 log(6V(27r) ). (13) 
As discovered by [C apacciol i, Caon fc D'Qnofriol ([19921 ) 
and I Graham fc Guzm an (200l), ETGs follow two differ- 
ent trends in the Re-(p)e plane, according to their lumi- 
nosity. The separation between the two families of bright 
and ordinary ellipticals occurs between Mb = — 19 and 
AIb = — 20. We adopt here a separation v alue of —20. By 
a linea r fit of the data in figure 9 of G raham fc Guzman! 
([2001 . we obtain a ^ 0.35 and /3 ^ —6.75 for the bright 
galaxies (Mb < -20) and a = -0.02 and /3 = 0.45 
for the ordinary ellipticals (Mb > —20). The latter 
val ue of a is consist ent with that of 0.34 ± 0.01 found 
bv lLa Barbera et al.l ([200 3b). who showed that the Ko- 
rmendy relation of bright ETGs does not change signif- 
icantly with redshift up to redshift z ~ 0.6 and that 
the intrisic scatter of the relation amount to 0.4 ± 0.03 
in (i.e. ~ 0.14 dex in R^^ g). In order to derive 
the Near-Infrared effective radius Re^, we use Eq. jT^j 
to compute Re^ g from Mb, and then transform R^j^ g 
into i?e£, . To this aim, we consider that ETGs have 
on average a radial color gradient of about —0.2 in 
B — K, and that their internal color gradients are ob- 
serv ed not to change significantly wi th galaxy l uminosity 
(see Peletier. Valentiin. fc Jamesonl [l 990: Pel etier et al.l 
I199C ). Following ISparks fc J5rgenseni (1993. ). the above 
value of the color gradient implies that the effective ra- 
dius of ETGs decreases by ~ 20% from B to K band. 
Thus, we derive Re^ from the relation 



Re 



O.SRe 



(14) 
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Fig. 1. — Physical consistency of the models. The logarithm of the distribution function / is plotted as a function of the relative 
binding energy £. The panels correspond to different values of the halo to stellar mass ratio, fi. From left to right and top to bottom, the 
panels correspond to = 0, 0.1, 1, 10, 10^, 10^. For each plot, as shown in the upper-left panel, curves with different colors correspond to 
different values of the Sersic index, while different line types denote different values of the ratio, x^, between the effective radii of the halo 



and stellar components. The /(e) has been computed by adopting natural units, where Mj^ 
was set to one. 

The Sersic parameter, n, of the stellar component 
depends on lum inosity through the magnitud e -Sersi c 
index relatio i i (ICaon. Capaccioli fc D'Onofriol [l993l) . 
iTruiillo et al.l ()2004[ ) presented this relation for a sam- 
ple of 200 ellipticals at redshift z ~ 0. A linear fit to the 
data in their figure 1 gives ^ 

logns = -0.1219- Ms - 1.6829, (15) 

where ub is the Sersic index of ETGs in the B-band. 
The Sersic index is not expected to change signifi- 
cantly from optical to NIR wavebands. For instance, 

^ We estimate the scatter of the luminosity— Sersic index re- 
lation from the dis tribution of points in Fig. 1 (right-panel) of 



1, Re 



1, and the gravitational constant 



ITruiillo et al.l l|2004l ). Assuming that, for a given magnitude, the 
smallest and largest Sersic index values mark the lower and upper 
2(T limits around the mean relation, we obtain a la dispersion of 
around 30% in ng at a given luminosity. 



as found by iLa Barbera et "aTl (|2008f ). ETGs have on av- 
erage \og{nr/nK) = — 0.007 ± 0.009, where and uk 
denote the r- and K-band Sersic indices. Hence, we set 
n = ub, and use Eq. [15] to derive also the NIR Sersic 
index of the stellar component. 

To express i^en as a function of Mb, we use the find- 
ing that dark matter halos follow a relation between the 
half-mass radius, Reo, and the average projected sur- 
face mass-density inside that radius, (ii)er,, similar to the 
Kormendy relation of galaxies (Graham et al. 2006; here- 
after GMM06). This resuh was obtained from GMM06 
for a sample of galaxy-sized dark-matter halos as massive 
as lO^^M©. We note that GMM06 derived the quantities 
i?eD and (/i)ei5 by fitting the projected halo density pro- 
filewith the Prugniel-Simien model (Prugniel & Simie3 
Il997( ). i.e. the same kind of profile as adopted here for 
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Fig. 2. — Stability of the models. The logarithm of the quantity g{r) (see Eq. [TTJ is plotted as a function of the logarithm of the 
dimensionless radius r/r^ i^. Colors and line types are the same as in Fig. [l] From left to right and top to bottom, the panels correspond 
to = 0, 0.1, 1, 10, 10^, 10^. Natural units have been adopted as for Fig. [T] 



the dark matter component of the models We write 

l0gi?e„ =<5- (Ai)e„ +C. (16) 

For systems more massive than 10^" Mq, GMM06 report 
a slope oi S ^ 1/3. This mass range corresponds to 
logi?ec > 0.4 (see fig. lb of GMM06). Performing a 
linear fit to the data in figure 2a of GMM06, we obtain 
C ~ 10/3, with i?ei3 being expressed in units of kpc. 

* We notice that GMM06 fitted the Prugniel-Simien model by 
treating the Sersic index as a free fitting parameter. Si nce we fix 
n = 3 for the dark-matter halo, the coefficients of Eq. 1161 taken 
from GMM06, might not be appropriate for our model calibration. 
When fitting a Sersic model with n = 3 to a Sersic profile with 
n = 3.6 (the average value found by GMM06), we find that the 
best-fitting eff'ective radius is ~ 20% smaller than the true value. 
However, due to the well-known correlation between effective ra- 
dius and mean surface brightness, this change in Re corresponds 
to a change in (i^e, such that points are moved a lmost parallel to 
the Kormendy relation l ILa Barbera et ani2003bl) . 



Then, we derive the mass of the dark-matter and stellar 
components as a function of the B-band magnitude, using 
the recent results obtained from lCappellari et al.l (|2006l ) 
(hereafter CAP06) for elliptical and lentic ular galaxies in 
the SAURON project (jBacon et al.ll200lD . From the re- 
lation between dynamical mass-to-light ratio in /—band 
and total mass of CAP06 (see their eq. 9), one obtains: 

M,^ + Afe„ = 1.175 • 10O-1317-O.528.A/« (^7) 

where M^^ and Mg^ denote the masses of the stellar 
and dark matter components within i?e^ . This rela- 
tion provides the total dynamical mass with an accuracy 
of ^ 30%. Following Fukugita, Shimasaku & Ichikaw^ 
(|1995[ ). we derive Eq. [17] by assuming a typical B — I 
color term ^ for elliptical galaxies of 2.23 and the B- and 

^ We notice that the assumption of a constant color 
term for ETGs is just a simplified assumption, since early- 
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I-band magnitudes of the Smi to be 5.51 and 4.08, re- 
spectively. Under the assumption of a radially constant 
M^/L ratio, one has = 2Me^. According to CAP06, 
Mf,j^ is about O.lGdex smaller ^ than the dynamical mass 
within Rej^ , i.e. Afg^^ 



Eq. [TTl one obtains: 



0.6918 (Afe^ + Me 



Thus, from 



Me^ = 0.81286 ■ 
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0. 1317-0. 528-Mi3 



(18) 



and 



Me, = 0.36214 •10"1317-0.528.M« 



(19) 

In order to relate Me^ to M, , we use the analytic ex- 
pression for the pro jected lumin osity profile of the Sersic 
model (see eq. 2 of lCiottilll999f) . Since the dark-matter 
component is described by a Sersic model having n — 3, 
we can write: 



Me, = . 7 6, 53 



Rp 



l/3^ 



(20) 



where 7 denotes the normalized incomplete gamma func- 
tion ^, and &3 — 5.6631. The quantity 63 is computed by 
setting n = 3 in the analytic approximation of b reported 
in Sec. 12.11 From Eqs. fTHland Eq. [TH one obtains: 



M„ = 



0.36214- 10 



0. 1317-0. 528-A/b 



1/3 



(21) 



In practice, for a given Mb, Re^ is computed from 
Eqs.fT2land[T4l and the quantities i?eo and M, are de- 
rived by solving simultaneouly Eqs. \n\ and [TBI This is 
equivalent to solve the non-linear equation 



7.1077H 



2.5(5 



-0.528MB+log 



7 6, 63 



Re 



Re 



1/3N 



55 
2. 



55 



(22) 

with respect to i^e, ■ We denote the first member of 



this equation as 9{Rea)- 
9{Reo) as a function of Re^, 



As an example. Fig. [3] plots 
for the case Mb = —21. 
The figure shows that Eq. [22] has in general two distinct 
solutions, corresponding to the points where the horizon- 
tal dashed line in the figure crosses the curve. One has a 



(and Me, < Me^ 



small-halo solution with iJe, < ReL 
and a large-halo case, whereby the dark-matter compo- 
nent is larger and more massive than the stellar one. In 
the small-halo case, the Me, value is four (eight) times 
smaller than Me^ for Mb = -22 (-20), while i?e, is 
three (ten) times smaller than _Rei, ■ This would imply 
that almost all the dark matter in ETGs should be en- 
closed within one i?ei , in disagr eement with dynamical , 
X-Ray, and weak lensing studies ([Matsushita et ai]ll998l : 

type systems ar e known to follo w a color-magnitude relation 
(e.g. [Visvanatha n &: Sandage|[T977l ). Though the above procedure 
can be generalized to account for a given color-magnitude relation, 
we decided to fix B — I. In fact, one should notice that the slope of 
the color-magnitune relation might be significantly affected from 
the aperture where color indices are d erived, due to th e existence 
of internal color gradients in galaxies l lScodeggiol I200l1) . with the 
slope flattening more and more as larger apertures are adopted. 

® This result was obtained under the assumption of a Kroupa 
IMF. 

The normalization of the incomplete gamma function is done 
by dividing it with the complete gamma function. We notice that 
in Sec. 2.1, we adopt a different notation where the 7 function is 
not normalized. 



[Wilson et all [20011: [Gerhard et al.ll2001[ ). Therefore, we 
consider here only the large-halo solutions of Eq. [22| We 
notice that Eq. [TCj applies to the case of massive galaxy- 
sized halos (Af, ^ 10^^), which might be appropriate 
only for bright galaxies {Mb < —20). For galaxies fainter 
than Mb = —20, we fix * the ratio of dark to stellar ef- 
fective radius to the value obtained for Mb = — 20 and 
then derive the total dark-matter mass from Eq. [211 

To summarize, we use the Kormendy and the 
luminosity-Sersic index relations to express R^j^ and n 
as a function of Mb- Then, by using Eg . [T8[ and solving 
Eq.[22l we also express M, , i?e, , and M^ as a function of 
AIb- In Tab.dl an example, we show the values of the 
five free parameters of the models that are obtained 
from the above procedure in six cases equally spanning 
the magnitude range of —22 to —17. In general, the 
procedure leads to have galaxy models where the dark 
matter component is less massive and less extended in 
lower luminosity systems. On the other hand, the rela- 
tive amount of dark matter within R^j^ does not depend 
on galaxy luminosity, in agreement with the finding of 
CAP06 (see Eqs. HI and [19] above). 

We remark that the above procedure derives the free 
parameters of the models by using the observed prop- 
erties of early- type systems at z ~ 0. Hence, one possible 
caveat when applying the above procedure to merging 
simulations is that such properties might not necessarly 
be the same for the high-redshift progenitors of ETGs. 
Moreover, one should consider that most of the observed 
relations (such as the Kormendy and the luminosity-size 
relations) of ETGs have significant intrinsic dispersion 
(see the values reported above), implying a dispersion, 
at a given magnitude, also in the parameter's values re- 
-^ported in Tab. [TJ 

l0g(i?e,) =0 

5. OPTIMAL SOFTENING LENGTH 

Performing discrete realizations of galaxy models re- 
quires that a given gravitational softening parameter, e, 
is adopted. The value of e should depend on the number 
of particles, N, defining the mass and spatial resolution 
of the simulation. Here, we discuss how to set e and N 
for the Sersic models. 

5.1. The optimal smoothing length 

Usually, the value of e is chosen with some ad hoc 
prescription. One fixes the total number of particles in 
the simulation (which is limited from the available CPU 
resources) and then assigns the e in or d er to achieve 
the desired spatial resolution. [Merritd ()1996( ) (here- 
after MER96) showed that the softening length of an 
N-body system can be chosen in an objective (optimum) 
way by minimizing the average error in the gravitational 
force computation over the whole space. Following a 
similar approach, we assign e by minimizing the aver- 
age error in the computation of the gravitational po- 
te ntial. We consider her e the spline softening kernel 
of iMonaghan fc Lattanziol (|T985), which is imple mented 
into the simulation code Gadget-2 (|Springel[[2005[ ) . Here- 
after, we express e in units of the effective radius, i?ei, ■ 

* Applying Eos. 1161 and 1221 also for Mb > —20 would lead to an 
improbable set of solutions where systems fainter than Mg = — 18 
would have dark-matter halos more massive than a galaxy with 
Mg = -22. 
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Re (kpc) 

Fig. 3. — Derivation of the e ffec tive radius of the dark matter component for a galaxy with Mb = —21. The i?eo is derived by solving 
the equation d{Re^ ) = (Eq. [22J. The horizontal dashed line marks the value of d{Ra^ ) = 0, while the vertical dashed line shows the 
effective radius iJe^ of the stellar component. The points of intersection between the horizontal line and the curve denote the values of 
iJeo which are consistent with our procedure. We consider only the large-halo solution (right part of the plot), with Rej^ > Rej^ (see the 
text) . 
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Mr. 
(IQiO Mq) 
(2) 


fief, 

(kpc) 
(3) 


M„ 

(lOW Mq) 

(4) 


Reo 

(kpc) 
(5) 


n 
(6) 


(IQlO Mq) 

(7) 


-22 


90.94 


14.84 


267.20 


108.60 


10.0 


20.25 


-21 


26.96 


5.07 


199.74 


75.50 


7.5 


6.00 


-20 


7.99 


1.73 


133.26 


45.50 


5.7 


1.78 


-19 


2.37 


0.90 


39.53 


23.79 


4.3 


0.53 


-18 


0.70 


0.87 


11.72 


22.81 


3.2 


0.16 


-17 


0.21 


0.83 


3.47 


21.88 


2.5 


0.05 



TABLE 1 

Derivation of the free parameters of the double Sersic models as a function of luminosity. The columns are: (1) B— band 

MAGNITUDE, Mg, OF THE STELLAR COMPONENT, (2) TOTAL STELLAR MASS, , (3) EFFECTIVE RADIUS OF THE STELLAR COMPONENT, Re^, 

(4) TOTAL MASS OF THE DARK MATTER HALO, M^ , (5) EFFECTIVE RADIUS OF THE DARK MATTER COMPONENT, fieo ■ (6) SeRSIC INDEX n, 

AND (7) MASS OF THE DARK MATTER HALO WITHIN Re^ ■ 



We start by considering the case of single Sersic mod- 
els. For a given Sersic index, n, and a given number of 
particles, N, we generate several realizations of the de- 
projected Sersic model. For a given realization, we calcu- 
late the softened gravitational potential at the position 
of each particle and the corresponding true gravitational 
potential (Eq. [5]). Then, the rms of the relative abso- 
lute differences between the softened and true potential, 
A(j>/(f), is computed over all the particles. We average the 
value of /S.4)/(f) over 100 realizations. Fig.|3]shows how the 
mean value of A<i>l (\) changes as a function of e. As exam- 



ple, the figure plots the case of a de Vaucouleurs model 
[n — 4) for two different values of iV. In both cases, there 
is a minimum in A(/)/0. Following an argument similar 
to that of MER96, the existence of a minimum can be 
explained as follows. For low e, the error is dominated by 
the differences between the point-like Newtonian poten- 
tial of each particle and the true gravitational potential. 
Increasing e, these differences become smaller and A^/(/) 
decreases. For large e, the discrete potential is smoothed 
on a scale larger than the typical interparticle separa- 
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e 



Fig. 4. — Mean value of the relative error on the gravitational 
potential as a function of the softenin length, e, for a Sersic model 
with n = 4. As shown in the upper-right corner of the plot, the 
dashed and solid curves correspond to discrete realizations with a 
different number of particles, N . 



tion ^ and the discrete potential is overly smoothed with 
respect to the true gravitational potential. Increasing e, 
this large-scale smoothing becomes more and more im- 
portant, and the value of A</)/0 increases as well. For 
a given number of particles, we define the position of 
the minimum as the optimal smoothing length, Eq. In- 
creasing the number of particles, the typical interparticle 
separation, djsi, decreases, and thus the optimal smooth- 
ing is obtained for smaller e. Fig. [5] plots as a function 
of N for different values of the Sersic index. The opti- 
mal softening length turns out to decrease as either n or 
N increase. This is due to the fact that, in both cases, 
the typical particle separation, dN, decreases. In par- 
ticular, when n increases, the mass profile of the model 
is more concentrated in the center and, at fixed iV, 
is smaller. As shown in Fig. [5l the trend of vs. N 
can be accurately modeled by a power law, eg — l3N~°', 
where both a and (3 depend on the value of n. The 
value of a changes from 0.28 for n = 1 to '--^ 0.54 for 
n = 7. For n < 2, the shape of the Sersic profile is 
flatter than for higher values of n, and the Eo is essen- 
tially proportional to the mean interparticle separation, 
with Eo cx iV^^/'^. For a de Vaucouleurs profile (n = 4), 
the value of a is ~ 0.4, in agreement with that of 0.44 
found by MER96 for the Hernquist model. For higher n, 
the Sersic profile becomes more and more peaked in the 
center and the value of a deviates more and more from 
the simple a ~ 1/3 expectation. Fig. [6] shows how the 
mean relative error on the potential, depends on 

the number of particles and the Sersic index when adopt- 
ing the optimal smoothing parameter. For a given Sersic 
model, the error decreases with N following the power- 
law A0/(/) cx N~^^'^, in agreement with what found by 
MER96 for the Hernquist model. For a given N, the 
error is larger for higher Sersic index. Hence, if a given 

^ The softening mostly affects the region where the potential 
changes more rapidly, i.e. the region inside the effective radius 
Rej^ ■ With typical interparticle separation, we refer to some sta- 
tistical estimator of the average particle-particle distance within 
that region, such as the mode or the median of the distribution of 
interparticle distances. 
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Fig. 5. — Dependence of the optimal softening length of one- 
component models, eg, on the number of particles, A'^, for different 
values of the Sersic index, n. Different colors correspond to differ- 
ent values of n as shown in the lower-left corner of the plot. Solid 
lines plot the best-fitted power laws to the trends of Co vs. N (see 
the text). The exponent a of each fitted power-law is reported on 
the top-right of the corresponding line. 
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Fig. 6. — The relative error on the gravitational potential, A(f>/<j), 
is plotted as a function of the number of particles, N, for three one- 
component Sersic models having n = 1 (dashed line), n = 4 (solid 
line) and n = 7 (dotted line), respectively. The gray line shows the 
power-law fit, A<l>/(j> oc N~^^'^ , to the points for n = 4. 

accuracy in the computation of the gravitational poten- 
tial has to be achieved, for higher n a larger number of 
particles has to be adopted. 

5.2. Models in isolation 

To perform discrete realizations of the models, one 
can adopt different softening lengths for the stellar and 
dark matter components, according to the optimal defi- 
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nition given above. However, these softening parameters 
represent an optimal choice only for one-component Ser- 
sic models, and we are not guaranteed that they provide 
also an accurate choice for the two-component models. 
To verify that the optimal prescription for e gives sen- 
sible results even in the case of two-component models, 
we compared the evolution of double and single Sersic 
models in isolation. As example, we consider here (1) 
a one-component model with Mj^ ~ 27 • 10^° Mq and 
i?ei, ~ 5 kpc, and (2) an S'^ model whose parameters 
are the same as those reported in Tab. [T] for the case 
Mb = —21. Model (1) is obtained by considering only 
the stellar component of model (2). To evolve the models 
in isolation, we adopt 50000 particles of luminous mat- 
ter in both cases and 75000 particles of dark matter for 
model (2). Looking at Fig. [6l we see that adopting the 
optimal smoothing parameter for these values of N al- 
lows an accuracy better than 10% on the gravitational 
potential to be achieved. The simulations were ran over 
5 Gyrs with the simulation code Gadget-2, using a Be- 
owulf system with thirty-two AMD-Opteron 244 proces- 
sors. As initial conditions, we created discrete realiza- 
tions of the models by computing their density profile 
and distribution function with the set of Fortran codes 
that are made publicly available (see App. [X)) . The soft- 
ening parameters were chosen according to Fig. [5l For 
the stellar component, we adopt eo = 0.013 kpc, while 
for the dark matter component we set Co = 0.053 kpc. 
Fig. [7] (upper panel) plots the relative absolute variation 
of the total energy of both systems, |A£'/£'o|, as a func- 
tion of time, where Eo is the total initial energy of the 
simulation. Apart from a small and slow secular drift, 
one can see that for both models the total energy of the 
system is preserved, with a value of \/S.E/Eq\ smaller 
than ~ 8% after 5 Gyrs. Fig. [7] (lower panel) also shows 
the evolution of the virial ratio, |2T/VF|, where T and W 
are the total kinetic and potential energy of the system, 
as a function of time. For both the single and S*^ models, 
the deviations from the virial equilibrium, 2T/W = 1, 
are small, amounting to at most ~ 0.7% in modulus af- 
ter 5 Gyrs. 

Fig. [8] plots, for the one-component model, the ra- 
dial profiles in mass, velocity dispersion, and anisotropy 
at T = Gyrs (left panels), and the relative variations 
of these profiles after the model has been evolved for 
5 Gyrs (right panels). The profiles are plotted in a ra- 



dial range of 



— 3eo to rraax = 5i?ei, • The value of 



''mm is chosen in order to avoid the inner region of the 
model which is affected by the smoothing in the gravita- 
tional potential. The maximum radius, r^ax, is set to a 
sensible value where one can compare the model to the 
observed profiles of ETGs. The simulation shows that 
the profile in mass is preserved within a few percentages 
over the whole radial extent. The velocity dispersion and 
the anisotropy profile are also preserved within ~ 10%. 
Fig. [9] plots the same profiles as in Fig. [8] for the stel- 
lar component of model (2). Remarkably, all the profiles 
are preserved even in this case within ~ 10% over at 
least 5 Gyrs. The same result was obtained when con- 
sidering the properties of the dark-matter component of 
model (2), and for all the 5*^ models whose parameters 
are listed in Tab. [1] 

6. SUMMARY AND DISCUSSION 
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Fig. 7. — Variations in total energy and virial ratio as a func- 
tion of time. Triangles and circles correspond to one- and two- 
component models, respectively (see the text). Notice that the 
deviations from conservation of total energy (|A_E/i?ol = 0) and 
the virial equilibrium (|2T/W| = 1) are small for both models. 
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Fig. 8. — Evolution in mass, velocity dispersion, and anisotropy 
profiles of single Sersic models. We plot the case of the model (1) 
described in the text. Left panels plot the mass (top), velocity dis- 
persion (middle) and anisotropy (bottom) profiles of the model at 
T = Gyr. The right panels show the relative absolute radial vari- 
ation of the profiles after T = 5 Gyrs. For each value of the spatial 
radius r, the variation is computed with respect to the initial value 
at that radius. We note that the variations from \2T/W\ = 1 are 
small for both models. 

We have presented models of ETGs consisting of a stel- 
lar component and a dark matter halo that follow the de- 
projected Sersic law. The models describe non-rotating, 
isotropic, spherical systems, whose density-potential pair 
is derived under the assumption that the stellar mass-to- 
light {Mj^/L) ratio of galaxies does not depend on radius. 

As mentioned in Sec. 12. 1[ the constant Mj^/L assump- 
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Fig. 9. — Same as Fig. |8] for the luminous component 
models. The case of model (2) is shown (see the text). 
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tion might not reflect the real physical properties of 
ETGs. Galaxies are observed to have internal color gra- 
dients, reflecting variations of stellar population proper- 
ties (such as age and m etallicity) from the galaxy cen- 
ter to the outskirts (e.g. iPeletier. Valentiin. fc JamesonI 
II99OD . It has been shown that (i) color gradi- 
ents are mainly driven by a mean metallicity gradi- 
ent in the range of = —0.2 to = —0.3, 
with an uncertainty of ~ 0.1; and that (ii) a small 
positive age gradient of Vt ~ 0.1 is also c onsis- 
tent with observations (see e.g. IPeletier et all [l990l: 
Saglia et"al] I2000t lldiart Michard fc de Freitas Pacheca 



20021: iLa Barbera et al.l l2003at iTamura fc Ohtal I2003D . 

Here, we denote as Vz and Vt the logarithmic variations 
of metallicity and age per decade in galaxy radius. From 
the theoretical viewpoint, age gradients are expected 
to arise in the formation of ETGs by gas-rich mergers, 
where early-type remnants are better described by a two- 
component stellar profile, with the two components hav- 
ing different ages (|Hopkins. Cox, fc Hernquistll2008f ) . We 
can use the above values of and V* to infer the corre- 
sponding radial variatio ns of JL. Using s i ngle s tellar 
populations models from'Bruzual fc CharlotI (|2003f ) with 
a Scalo IMF and an age of 12Gyr one obtains that a 
metallicity gradient of Vz = —0.2 (—0.3) corresponds to 
a variation of 34% (51%) in the B-band / L per decade 
of galaxy radius. This variation largely decreases in K- 
band, where the inferred variation of M^/ L amounts to 
~ 15% (24%). Considering a positive age gradient of 
O.ldex, the M^jL variation would further decrease to 
about 7% (10%) in K-band, while the above uncertainty 
on color gradients would translate to an error of about 
one third in the estimated M^jL percentages. We con- 
clude that, provided one adopts the K-band light profile 
of ETGs to infer the underlying distribution of stellar 
matter, the assumption of a constant M^/L is empiri- 

In a cosmology with f2m = 0.3, Ca = 0.7, and Ho = 70 km 
s~ ^ Mpc~ ^ , this would correspond to a formation redshift of z ~ 4. 



cally well motivated. 

For what concerns the other assumptions underlying 
the S"^ models, one should notice that ETGs actually 
span a wider range of kinematical and structural prop- 
erties than that considered here. For instance, the 
models populate the origin of the anisotropy {v/a vs. el- 
lipticity) diagram, while ETGs populate different regions 
of it. In order to explore the corresponding effect on dry- 
merging simulations, some studies have realized merg- 
ing simulations wher e the p rogenitors are obtained by 
either dissipatio nless (|Naab. Khochfar, fc Burkert 2006[) 
or dissipational (|Cox et al.ll2006rrRobertson et aLlbood ) 
merging of disk systems. This re-merger approach 
has the main advantage that progenitors span a wide 
range of ETC properties, such as w/cr, ellipticty, and 
isophotal shape. Though neglecting these aspects, 
the S*^ models have the main advantage of allow- 
ing one to explore a key observational feature: the 
wide r ange of profile shapes observed in ea rly-type sys- 
tems ("Ca on. Capaccioli fc D'Onofriol n"993f ). Moreover, 
re-merging of 5'^ would likely allow one to further en- 
large the range of kinematic and isophotal properties of 
merging progenitors. 

The free parameters of the two components of S"^ mod- 
els are assigned in order to match the observed properties 
of ETGs as well as recent results of N-body simulations 
of galaxy-sized dark matter halos. We report a concise 
reference to the basic integral equations that define the 
density-potential pair and the distribution function of 
the deprojected Sersic law, showing how these equations 
can be used to define the S*^ models. We show that for 
all possible values of the free parameters of the models, 
the total distribution function is always non-negative de- 
fined, implying that the models are physically admissible 
solutions of the coUisionless Boltzmann equation. More- 
over, the first derivative of the total distribution function 
is always non-negative defined, implying that the models 
are stable against radial and non-radial perturbations. 
For a given Sersic model, we present an objective pre- 
scription to adopt an optimal smoothing length of dis- 
crete model realizations. The optimal smoothing length 
is defined as the softening parameter that minimizes the 
error on the gravitational potential of the system, and 
depends on the Sersic index n as well as on the number 
of particles of the simulation. The power-law relations 
that describe these trends are reported, with the aim of 
providing a prescription to create discrete realizations of 
S*^ systems, whose discrete gravitational potential closely 
matches the true model potential. As a caveat, when 
using such a prescription for merging simulations, one 
should notice that the optimal smoothing length for the 
progenitors might not necessarely concide with the op- 
timal softening for the merging remnants, depending on 
the structural properties (i.e. the Sersic index) of the 
merging end-products. This issue can be addresses by 
exploring the effect of changing the number of particle, 
and the corresponding smoothing length, of the colliding 
systems. 

We provide the Fortran code that allows one to calcu- 
late all the properties of single and double Sersic models. 
The code together with the recipes for computing the op- 
timal softening scale are intended as general tools to per- 
form merging simulations of early-type galaxies, whereby 
the structural non-homology of these systems (i.e. the 
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variation of the shape parameter along the galaxy se- 
quence) might be taken into account. In a compan- 
ion contribution (Coppola et al. 2009b, in preparation), 
we use the models to investigate how dissipation-less 
(major and minor) mergers affect the structural proper- 
ties of ETGs, such as the shape of their light profile and 
their stellar population gradients. 



We thank L. Mayer and E. D'Onghia for the help- 
ful comments and suggestions. We also thank the ref- 
eree who provided several comments/suggestions which 
helped us to significantly improve this manuscript. 



APPENDIX 
FORTRAN CODES 

The properties of both the single and double Sersic models are computed by a set of FORTRAN routines. All the 
Fortran codes are made publicly available For the one-component models, the code allows the user to calculate 
the density, mass, and gravitational potential profiles (by a numerical integration of Eqs. O 21 and [5]), as well as the 
distribution function (App. [B|. Other quantities, such as the total potential and gravitational energy of the system, 
its spatial and projected velocity dispersion profiles, are also computed by specific Fortran routines. For the double 
Sersic model, since the computation of the density-potential pair is time-demanding, we proceed as follows. 

- For a given value of the Sersic index n, that characterizes the luminous component of the model, we calculate 
the dimension-less mass, density, potential and the first and second derivatives of the density profile over a grid 
in the dimension-less spatial radius x. The same computation is done for the Sersic index of the dark matter 
component, n = 3 (Sec. 12. 3p . 



- The total density-potential pair and the distribution function are then obtained by interpolating the above radial 
profiles. To this effect, the values of the parameters /i and xd of the model have to be provided (Sec. 12. 3p . 

The software to perform this interpolation procedure is also provided. 

DISTRIBUTION FUNCTION OF THE DOUBLE SERSIC MODEL 

,2 

In order to apply the Eddington inversion (Eq. [TU)) . one has to calculate the function -j^, where p is the spatial 
density profile and ^E* = —ip + ipo is the rescaled gravitational potential (see Sec. 12. 3p . We start from the following 
identity: 



fp_ 



(Pp 



Then, using the fact that p{r) = p^. 



dr J 

and if{r) = ip^ 



dp 

dr 



dr J dr^ 



(Bl) 



ip^ , one obtains the following expression: 



■fp 


m 


-m 


SI 








dr'^ 






dr^ 




dr^ dr 


dr dr'^ 



d^PL dVr. _ dp^ d'^Lp^ 
dr^ dr dr dr^ 



dPr^ d'^iPo _ dp^ 
dr dr^ dr dr"^ 



(B2) 



The first and second derivatives of p^ and p^ can be derived by numerically differentiating Eqs.[3]and[6] The derivatives 
of the gravitational potential and the density profile can be obtained from the expression of the gravitational potential 
and the mass profile of the stellar and dark matter components, using the following identities: 



d^L 
dr 



2 \x=r/R,,^ ) 



dp^ _ GM^ p M{x/x^) 



dr^ 



G 



dr 

= AttG 
p 



M, 



2 [x/x^Y 
62" 



Rl^ 27rnr(2n 



i?3 X^ 



dPr, 

dr 

dPo ^ 
dr 



mT{2m) 
M, fo2« 



ipD{x/x^) 
dp . 



2 Id^, 
X dr ' 

2M(x/xJ 
{x/x^Y 



Ri^ 27mr(2n) dx '^"''/^^ 

p 62m ^~ 



xl 27rmr(2m) dx 



(B3) 
(B4) 
(B5) 
(B6) 
(B7) 
(B8) 
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u2n 



dr'^ Rl 27rnr(2n) dx= 



and 



dr^ R^^ 27rTOr(2m) dx 



2 \ x—x / X jj : 



(B9) 



(BIO) 



These equations show that the f{£) is completely defined by the first and second derivatives of the density profile, 
the gravitational potential and the mass profiles of the two Sersic components. In order to calculate f{£), we derive 

numerically the functions p, 0, and M, and then, using Eg. IB2( we evaluate Eg. 1101 

To prove the stability of the double Sersic models, one has to prove the condition ^ > (see Sec. [3]). From the 

Eddington formula, a necessary condition is 4^ > 0. From Eg. IBll this condition can be written as 



fp 
dr"^ 



dr 



dp 
dr 



d-^ 
dr 



dr^ 



d* 

dr 



fp 
dr"^ 



d* 
dr 



d^^ 
dr^ 



> 



(Bll) 



Since ^ is negative (i.e. the gravitational potential is a monotonically increasing function of r), the previous condition 



is eguivalent to 



as stated in Sec. 



g{r;n,n,x^) 



dr^ 



dr 



dr^ 



>0, 



(B12) 
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